Imu calibration

ABSTRACT

A method of calibrating an inertial measurement unit, the method comprising: (a) collecting data from the inertial measurement unit while stationary as a first step; (b) collecting data from the inertial measurement unit while repositioning the inertial measurement unit around three orthogonal axes of the inertial measurement unit as a second step; (c) calibrating a plurality of gyroscopes using the data collected during the first step and the second step; (d) calibrating a plurality of magnetometers using the data collected during the first step and the second step; (e) calibrating a plurality of accelerometers using the data collected during the first step and the second step; (f) where calibrating the plurality of magnetometers includes extracting parameters for distortion detection and using the extracted parameters to determine if magnetic distortion is present within a local field of the inertial measurement unit.

CROSS REFERENCE TO RELATED APPLICATIONS

The present application claims the benefit of U.S. Provisional PatentApplication Ser. No. 62/268,175, titled “IMU CALIBRATION,” filed Dec.16, 2015, the disclosure of which is incorporated herein by reference.

INTRODUCTION TO THE INVENTION

The present disclosure is directed to methods and techniques forcalibrating inertial measurement units (IMUs) and, more specifically, tomethods and techniques that involve calibrating IMUs to account forlocal environments within which the IMUs operate.

It should be noted that Patent Cooperation Treaty applicationPCT/US14/69411, filed Dec. 9, 2014 is hereby incorporated by referenceand appended hereto as Appendix A.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a graphical depiction of a ferromagnetic device mounted to aninertial measurement unit and being rotated around an X-axis.

FIG. 2 is a graphical depiction of the ferromagnetic device mounted tothe inertial measurement unit of FIG. 1, being rotated around a Y-axis.

FIG. 3 is a graphical depiction of the ferromagnetic device mounted tothe inertial measurement unit of FIG. 1, being rotated around a Z-axis.

DETAILED DESCRIPTION

The exemplary embodiments of the present disclosure are described andillustrated below to encompass inertial measurement units and methods ofcalibrating the same. Of course, it will be apparent to those ofordinary skill in the art that the embodiments discussed below areexemplary in nature and may be reconfigured without departing from thescope and spirit of the present invention. However, for clarity andprecision, the exemplary embodiments as discussed below may includeoptional steps, methods, and features that one of ordinary skill shouldrecognize as not being a requisite to fall within the scope of thepresent invention.

Many current devices make use of integrated sensor components to captureand track information related to the use, location, and status of thedevice. One such system, the inertial measurement unit (IMU), is used totrack object motion. A common packaging of an IMU comprises a singleintegrated circuit (IC) consisting of multiple inertial sensors, whichcan be a combination of accelerometers, gyroscopes, and magnetometers.These three sensor types can also be separated into individual ICspopulated on the same circuit board. These two cases (single IC andmultiple ICs) can be considered the same for purposes of the instantdisclosure.

For all such sensors of the IMU, a calibration step may be performed toensure the proper referencing and operation in the local environment.For example, accelerometers may be calibrated to determine a referencewith respect to the direction of gravity and gyroscopes may becalibrated to some known rotational velocity. Magnetometer calibrationis more complex, as the process requires knowledge of the local magneticfield strength and direction at a multitude of magnetometerorientations. The local environment calibration can be carried out byhand or by machine. The minimum considerations for local environmentcalibration should involve maneuvering of the IMU around the threeorthogonal axes of the IMU, and a stationary period.

For applications requiring highly reliable and reproducible calibrationof an IMU, a calibration machine may be used. A calibration machine maybe constructed so that it does not introduce additional magneticdistortion to the IMUs. The calibration machine should produce movementfor calibration for all the sensors of the IMU as well as movement inorder for an algorithm to inspect the quality of the calibration.

The process for calibrating measurements from uncalibratedaccelerometers, gyroscopes, and magnetometers of an IMU is describedhereafter. This exemplary process is applicable to a single IMU ormultiple IMUs populated on the same circuit board, where all sensors(uncalibrated accelerometers, gyroscopes, and magnetometers) may becalibrated simultaneously.

In exemplary form, the calibration procedure may include the followingroutines to allow adequate data to calculate and verify the calibration.This exemplary calibration comprises four stages that include, withoutlimitation: (1) calibrating the gyroscopes; (2) calibrating themagnetometers routine; (3) calibrating the accelerometer routine; and,(4) calibrating the checking/validation routine.

Calibrating the gyroscopes of each IMU should provide for the gyroscopesbeing stationary for a set amount of time. As part of calibrating thegyroscopes, the instant process is determining the bias of the gyroscopeby arithmetic mean as set forth in Equation 1 immediately below:

$\begin{matrix}{{{{Gyro}_{{mean}_{i}}\lbrack\ldots\rbrack} = {\overset{\_}{{DATA}_{{gyro}_{\iota}}}}_{{gyro}\mspace{14mu} {calibration}}},{i = 1},2,\ldots} & {{Equation}\mspace{14mu} 1}\end{matrix}$

Calibrating the magnetometers of each IMU should involve rotating theIMU around all of its orthogonal axes. In this fashion, thismagnetometer calibration process first determines the soft and hard irondistortion using multi-dimensional Ellipsoidal Fitting function (EF_fun)set forth in Equation 2 immediately below:

$\begin{matrix}{{SFe}_{{orig}_{i}},{{HFe}_{{orig}_{i}} = {EF}_{{fun}{({{DATA}_{{mag}_{i}}{\lbrack{{{Mag}\mspace{14mu} {calibration}\mspace{14mu} {being}}->{end}}\rbrack}})}}},\mspace{20mu} {i = 1},2,\ldots} & {{Equation}\mspace{14mu} 2}\end{matrix}$

where: SFe_(orig) is the transformation of ellipsoid to sphere, andHFe_(orig) is the sensor and local bias.

The resultant data from Equation 2 are then processed using Equation 3,shown immediately below:

$\begin{matrix}{{procMAG}_{{orig}_{i}} = {{processMag}( {{{DATA}_{{mag}_{i}} \quad{\lbrack {{{Acc}\mspace{14mu} {calibration}\mspace{14mu} {begin}}->{{end}\mspace{14mu} {of}\mspace{14mu} {checking}\mspace{14mu} {routine}}} \rbrack,{SFe}_{{orig}_{i}},{HFe}_{{orig}_{i}}} )},{i = 1},2,\ldots}\mspace{14mu} )}} & {{Equation}\mspace{14mu} 3}\end{matrix}$

Thereafter, the alignment differences among magnetometers arecalibrated. In doing so, Equation 4 may be utilized, with the assumptionthat all magnetometers are aligned to a reference alignment using thisequation. Alternatively, the magnetometers may be aligned to a differentreference alignment if desired. In Equation 4, TMAG_(ij) represents thealignment transformation of magnetometer “j” to reference magnetometer“i”.

$\begin{matrix}{{{TMAG}_{ij} = {{procMAG}_{{orig}_{i}}\backslash {procMAG}_{{orig}_{i}}}},{i = 1},{j = 2},{3\mspace{14mu} \ldots}} & {{Equation}\mspace{14mu} 4}\end{matrix}$

Calibrating the accelerometers of each IMU should involve rotating theIMU around two axes perpendicular to gravity. In this fashion, theaccelerometer calibration is calculated with the multi-dimensionalEllipsoidal Fitting function (EF_fun) set forth immediately below asEquation 5:

$\begin{matrix}{{AccScale}_{{orig}_{i}},{{AccBias}_{{orig}_{i}} = {EF}_{{fun}{({{DATA}_{{Acc}_{i}}{\lbrack{{{Acc}\mspace{14mu} {calibration}\mspace{14mu} {begin}}->{end}}\rbrack}})}}},{i = 1},2,\ldots} & {{Equation}\mspace{14mu} 5}\end{matrix}$

Calibrating the checking routine should involve rotating the IMU aroundthe axis collinear with gravity. The result is best if the motion of thechecking routine is performed multiple times at different angles withrespect to gravity.

After calibrating the sensors, one may take into account manufacturingand IC assembly placement errors to extract or account for these errors.In this fashion, the instant disclosure includes a calibration sequencetaking into account any mal-alignment among sensors, which arecalibrated according to Equation 6-10 (applying the calibrationparameters) provided immediately below. In doing so, all sensor datacollected during the calibration may be subject to these calibrationequations/parameters.

procGyro_(i)=processGyro(DATA_(gyro) _(i) ,Gyro_(bias) _(i) ),i=1,2,  Equation 6

procAcc_(i)=processAcc(DATA_(Acc) _(i) ,ACC _(bias) _(i) ),i=1,2,  Equation 7

procMAG _(i)=processMag(DATA_(mag) _(i) ,SFe _(orig) _(i) ,HFe _(orig)_(i) ), i=1,2,  Equation 8

Continuing the example as aligning all magnetometer to the first one,

alignProcMAG _(j)=procMAG _(j) *TMAG _(ij) , i=1, j=2,3  Equation 9

alignProcMAG_(i)=procMAG _(i)  Equation 10

With multiple IMUs rigidly fixed together, the aligned and calibratedsignals should be roughly equal among each sensor type in the case of nonoise or isotropic (to the sensor) noise.

Despite calibration with the local environment, magnetometers areespecially prone to transient local disturbance. Specifically,introduction of additional ferromagnetic sources or movement offerromagnetic sources after initial calibration (described previously)can distort the sensor signal and result in inaccurate readings.Distortion states can be classified as hard iron distortion, which isprimary caused by the introduction of ferromagnetic materials such assteel tools that can deflect the field lines and produce an offset onthe magnetometer sensors, and soft iron distortion, which is caused byan irregular magnetic field (such as magnets) that can alter the shapeof the magnetic field. Given that these devices causing distortions willenter and exit the application field (area of use) of the IMU atunpredictable orientations, positions, and times, it is desirable thatthese distortions be detected and accounted for. By way of example, onearea where distortions are commonplace is an operating room during asurgical procedure. More specifically, surgical instruments within anoperating room are frequently repositioned within the operating field atdifferent times and locations depending on the stage of the surgicalprocedure, which necessarily causes magnetic field distortions. Itshould be noted, however, that the exemplary disclosure is not limitedto calibration of IMUs for surgical navigation or any particular fieldof use. Rather, the instant disclosure is applicable to calibrating anyIMU in any field of use where local distortions may be present.

There are two essential parts to detect transient magnetic distortion.The first part is that the IMU must comprise two or more magnetometers(that may be populated on a common circuit board). The magnetometers arecalibrated and aligned based the foregoing calibration teachings. Thesecond part is the detection algorithm, which consists of two processes:(1) extracting parameters for distortion detection, which are calculatedfrom the calibration data; and, (2) use of extracted parameters oncaptured IMU data to check if distortion is present.

What follows is a detailed explanation for calculating appropriatedistortion threshold parameters. In general, one may use an algorithmfor this calculation that may make use of the following assumptions: (1)IMU motion (or lack of motion) is detectable by all IMUs (accelerometer,gyroscope and magnetometer); (2) the angle and length of the calibratedmagnetometers should have approximately identical values; (3) thecalibrated magnetometer vectors should be approximately unity inmagnitude; (4) the angle between the calibrated magnetometer vector andcalibrated accelerometer vector should not change; (5) the radius of thecalibrated magnetometers should be approximately equal to one; and, (6)the quaternions calculated using different magnetometers should beapproximately equal.

Sufficient deviation from assumption 1 is deemed a motion anomaly anddeviations from assumptions 2-6 are deemed a magnetic anomaly.Parameters for motion detection are calculated for each IMU sensor(accelerometer, gyroscope, and magnetometer) using data collected duringlocal environment calibration. Using stationary and dynamic data,appropriate thresholds of motion are calculated for every IMU sensor, sothat the calibrated IMU sensor output above the threshold indicates thatthe sensor is in motion and below which the sensor can be consideredstationary (see Equation 2). A motion anomaly is considered to occur asin Table (2), which is described in more detail hereafter. In short, ifa magnetometer has detected motion, but other sensors do not, a motionanomaly is said to exist.

For the magnetic anomaly, a deviation from assumption 2 is used tocreate a cost function, which is calculated as the magnitude of thevector difference weighted with the angular derivation between thevectors (see Algorithm 1 and Algorithm 4). During local environmentcalibration, expected values for the deviation from unity and fromidentical signals are calculated to form a cost function. The thresholdis calculated from the value of this cost function during localenvironment calibration (when no additional distortion is presented).Then, as part of an anomaly detection process, any values of this costfunction above the threshold are considered a magnetic anomaly.

For assumption 3, the magnetic strength is calculated as the mean ofcombined magnitudes (Algorithm (5)), which without noise is expected tobe approximately 1. Then, as part of the anomaly detection process, anyvalue outside the calculated magnetic strength (including a presettolerance) is considered an anomaly.

For assumption 4, the derivation is calculated as the angular derivationbetween the magnetometer and the accelerometer vectors (Algorithm (6)).Then, as part of the anomaly detection process, any value outside thecalculated angle (including a preset tolerance) is considered ananomaly.

For assumption 5, the deviation is calculated by using all of themagnetometer readouts, along with a preset of artificial points, toestimate a radius of the magnetic field using an ellipsoid estimationfunction as shown in Algorithm (7). The artificial points are used tosupport the estimation calculation. However, if the IMU contains morethan 11 magnetometers, no artificial point is needed to resolve theradius. Then, as part of the anomaly detection process, any valueoutside the calculated radii (including a preset tolerance) isconsidered an anomaly.

For assumption 6, the deviation is determined by first calculating theorientation of the IMU using sensor fusion algorithm (e.g. ExtendedKalman Filter) using different magnetometers as input, and an angulardifference among the output orientation is calculated as shown inAlgorithm (8). Then, as part of the anomaly detection process, anyvalues greater than the angular different with a preset tolerance isconsidered an anomaly.

If motion or magnetic anomaly occurs, some unaccounted for magneticdistortion is present. An alert system for distortion provides feedbackindicating the reliability of the orientation reported by the sensorsmay be inaccurate. What follows is a more detailed discussion of thealgorithms utilized to detect motion and magnetic anomalies.

The anomaly detection is based on three categories of detection. Thefirst is sensor fault detection, the second is magnetometer radiiverification, and the third is quaternion output discrepancies. Sensorfault detection consists of the following functions: (i) motion andmagnetic anomalies; (ii) magnetic strength; and (iii) angle betweenmagnetometer and accelerometer vectors. The following method is used toconsistently extract threshold values (TVs) for motion anomalies andcost functions from the calibration. The method utilizes Algorithm (1)as follows:

  Algorithm  (1)Iteration:  for  [beginning  of  gyro  calibration → end  of  checking  routines]$\mspace{79mu} {{{{Gyro}_{{mean}_{i}}\lbrack\; \ldots \;\rbrack} = {{\overset{\_}{procGyro}}_{{current} - {kernel}}^{current}}},{i = 1},2,\ldots}$Gyro_(std_(i))[ … ] = std(procGyro[current − kernel → current]), i = 1, 2, …Acc_(std_(i))[ … ] = std(procAcc_(i)[current − kernel → current]), i = 1, 2, …Mag_(std_(i))[ … ] = std(alignProcMAG_(i)[current − kernel → current]), i = 1, 2, …$\mspace{79mu} {{{{lenDiff}_{i,j}\lbrack\; \ldots \;\rbrack} = \frac{\sqrt{\sum\limits_{axis}( {\sum\limits_{{current} - {kernel}}^{current}\begin{pmatrix}{{alignProcMAG}_{i} -} \\{alignProcMAG}_{j}\end{pmatrix}^{2}} )}}{kernel}},\mspace{20mu} {i \neq j},{i = 1},{2\mspace{14mu} \ldots}\mspace{14mu},{j = 1},{2\mspace{14mu} \ldots}}$$\mspace{20mu} {{{{angDiff}_{i,j}\lbrack\; \ldots \;\rbrack} = \frac{{1 - \begin{pmatrix}{{\overset{\_}{alignProcMAG}}_{l_{{current} - {kernel}}}^{current} \cdot} \\{\overset{\_}{alignProcMAG}}_{J_{{current} - {kernel}}}^{current}\end{pmatrix}}}{\begin{matrix}{{{norm}( {\overset{\_}{alignProcMAGcur}}_{l_{{current} - {kernel}}}^{current} )} \cdot} \\{{norm}( {\overset{\_}{alignProcMAG}}_{J_{{current} - {kernel}}}^{current} )}\end{matrix}}},\mspace{20mu} {i \neq j},{i = 1},{2\mspace{14mu} \ldots}\mspace{14mu},{j = 1},{2\mspace{14mu} \ldots \mspace{20mu} {{{metricLenAng}_{i,j}\lbrack\; \ldots \;\rbrack} = {{{lenDiff}_{i,j}\lbrack\; \ldots \;\rbrack}*{{angDiff}_{i,j}\lbrack\mspace{11mu} \ldots \mspace{11mu}\rbrack}}}},\mspace{20mu} {i \neq j},{i = 1},{2\mspace{14mu} \ldots}\mspace{14mu},{j = 1},{2\mspace{14mu} \ldots {{{currentMagVal}\;\lbrack\; \ldots \;\rbrack} = \frac{1}{{\max (i)}( {\sum\sqrt{\begin{matrix}{\sum\limits_{{current} - {kernel}}^{current}\;} \\( {{alignProcMAG}_{l}*{alignProcMAG}_{J}} )\end{matrix}}} )}}},\mspace{20mu} {i \neq j},{i = 1},{2\mspace{14mu} \ldots}\mspace{14mu},{j = 1},{2\mspace{14mu} \ldots}}$  End

Once the iteration from Algorithm (1) is completed, the mu and sigma ofthe cost function can be determined using Equations 11 and 12,immediately recited below.

μ=currentMagVal _(mag and acc calibration)  Equation 11

δ=SSF1*std(currentMagVal)_(mag and acc calibration)  Equation 12

Where: (a) SSF1 is a sensitivity scale factor; and,

(b) δ serves as the detection tolerance based on standard derivation(i.e., how many standard deviations away from the mean is consideredacceptable).

Algorithm (2) is utilized to calculate motion and no motion (stationary)states for the sensors.

$\begin{matrix}{{{Stationary}_{{gyro}_{std}} = {\overset{\_}{{Gyro}_{std}}}_{{gyro}\mspace{14mu} {calibration}}}{{Motion}_{{gyro}_{std}} = {\overset{\_}{{Gyro}_{std}}}_{{acc}\mspace{14mu} {calibration}}}{{Stationary}_{{gyro}_{mean}} = {\overset{\_}{{Gyro}_{mean}}}_{{gyro}\mspace{14mu} {calibration}}}{{Motion}_{{gyro}_{mean}} = {\overset{\_}{{Gyro}_{mean}}}_{\; {{acc}\mspace{14mu} {calibration}}}}{{{Stationary}_{{Acc}_{i}} = {\overset{\_}{{Acc}_{\iota_{std}}}}_{{gyro}\mspace{14mu} {calibration}}},{i = 1},2,\ldots}{{{Motion}_{{Acc}_{i}} = {\overset{\_}{{Acc}_{\iota_{std}}}}_{{acc}\mspace{14mu} {calibration}}},{i = 1},2,\ldots}{{{Stationary}_{{MAG}_{i}} = {\overset{\_}{{MAG}_{\iota_{std}}}}_{{gyro}\mspace{14mu} {calibration}}},{i = 1},2,\ldots}{{{Motion}_{{MAG}_{i}} = {\overset{\_}{{MAG}_{\iota_{std}}}}_{{Checking}\mspace{14mu} {routine}}},{i = 1},2,\ldots}} & {{Algorithm}\mspace{14mu} (2)}\end{matrix}$

A check is performed here to ensure all motion states should have valuesgreater than the corresponding stationary states.

Algorithm (3) is utilized to calculate TVs for motion anomalies.

$\begin{matrix}{{{{TV}_{{gyro}_{std}} = {{{SSF}\; 2*( {{Motion}_{{gyro}_{std}} - {Stationary}_{{gyro}_{std}}} )} + {stationary}_{{gyro}_{std}}}}{{TV}_{{gyro}_{mean}} = {{{SSF}\; 3*( {{Motion}_{{gyro}_{mean}} - {Stationary}_{{gyro}_{std}}} )} + {Stationary}_{{gyro}_{mean}}}}{{{TV}_{{Acc}_{i}} = {{{SSF}\; 4*( {{Motion}_{{Acc}_{i}} - {Stationary}_{{Acc}_{i}}} )} + {Stationary}_{{Acc}_{i}}}},{i = 1},2,\ldots}{{{TV}_{{MAG}_{i}} = {{{SSF}\; 5*( {{Motion}_{{MAG}_{i}} - {Stationary}_{{MAG}_{i}}} )} + {Stationary}_{{MAG}_{i}}}},{i = 1},2,\ldots}\mspace{20mu} {SSF}\; 2\mspace{14mu} {to}\mspace{14mu} {SSF}\; 5\mspace{14mu} {are}\mspace{14mu} {the}\mspace{14mu} {sensitivity}\mspace{14mu} {scale}\mspace{14mu} {{factor}.{It}}\mspace{14mu} {is}\mspace{14mu} {typically}\mspace{14mu} {set}\mspace{14mu} {between}\mspace{14mu} 0\mspace{14mu} {and}\mspace{14mu} 1},{{and}\mspace{14mu} {should}\mspace{14mu} {be}\mspace{14mu} {adjusted}\mspace{14mu} {based}\mspace{14mu} {on}\mspace{14mu} {sensors}\mspace{14mu} {and}\mspace{14mu} {desired}\mspace{14mu} {sensitivity}\mspace{14mu} {of}\mspace{14mu} {the}\mspace{14mu} {detection}\mspace{14mu} {{method}.}}} & {{Algorithm}\mspace{14mu} (3)}\end{matrix}$

Determination of a magnetic anomaly may be based upon a cost functionTV, which is accomplished by calculating the cost function over theperiod of the calibration period as set forth in Algorithm (4).

  Algorithm  (4)Iteration:  for  [beginning  of  gyro  calibration → beginning  of  Check  1]${{pdfVal}\;\lbrack\; \ldots \;\rbrack} = {{\frac{1}{\sqrt{2\; \pi}*\delta}e^{{- 0.5}*{(\frac{({1 - \mu})}{\delta})}^{2}}} - {\frac{1}{\sqrt{{2\; \pi}\;}*\delta}e^{{- 0.5}*{(\frac{({{{currentMagVal}\;\lbrack\; \ldots \;\rbrack} - \mu})}{\delta})}^{2}}}}$${{{costFunction}\;\lbrack\; \ldots \;\rbrack} = {S\; S\; F\; 6*{{{{pdfVal}\;\lbrack\; \ldots \;\rbrack}*\frac{\sum{{metricLenAng}_{i,j}\lbrack\mspace{11mu} \ldots \mspace{11mu}\rbrack}}{\max (i)}}}}},\mspace{20mu} {i \neq j},{i = 1},{2\mspace{14mu} \ldots}\mspace{14mu},{j = 1},{2\mspace{14mu} \ldots}$  endTV _(costFunction)=max(costFunction)+(costFunction)+std(costFunction)

SSF6 is a scale factor that is used to prevent losing the precision ofthe cost function from the small number.

In order to consistently extract TVs for discerning the strength of amagnetic anomaly from the calibration data, one can use Algorithm (5).

  Algorithm  (5)Iteration:  for  [beginning  of  mag  calibration → end  of  mag  calibration]$\mspace{20mu} {{{{MatStr}_{i}\lbrack\; \ldots \;\rbrack} = {{\sqrt{\sum\limits_{axis}( {alignProcMAG}_{i} )}\mspace{14mu} i} = 1}},2,\ldots}$  end   min_MagStr_(i) = min (MagStr_(i))  i = 1, 2, …  max_MagStr_(i) = max (MagStr_(i))  i = 1, 2, …

In order to consistently extract TVs for discerning the angle betweenmagnetometer and accelerometer vectors from the calibration data, onecan use Algorithm (6).

  Algorithm  (6)Iteration:  for  [beginning  of  gyro  calibration → end  of  checking  routine]  Use  Extended  Kalman  Filter  to  estimation  the  orientation  during  the  calibration  q_(i)[ … ] = EKF(procGyro, procAcc, alignProcMAG_(i)),   i = 1, 2, …$\mspace{20mu} {{{e_{i}\lbrack\; \ldots \;\rbrack} = \begin{bmatrix}{q\; 2_{i}} & {q\; 3_{i}} & {q\; 4_{i}}\end{bmatrix}},{i = 1},2,\ldots}$$\mspace{20mu} {{{{\hat{e}}_{i}\lbrack\; \ldots \;\rbrack} = \begin{bmatrix}0 & {{- q}\; 4_{i}} & {q\; 3_{i}} \\{q\; 4_{i}} & 0 & {{- q}\; 2_{i}} \\{{- q}\; 3_{i}} & {q\; 2_{i}} & 0\end{bmatrix}},{i = 1},2,\ldots}$  Rp_(i)[ … ] = (q 1_(i)² − e_(i)^(′)e_(i)) * I_(3 × 3) + 2(e_(i)e_(i)^(′)) − 2q 1_(i)ê_(i), i = 1, 2, …  MF_(i)[ … ] = Rp_(i) * alignProcMAG_(i), 1 = 1, 2, …  GF_(i)[ … ] = Rp_(i) * procAcc, i = 1, 2, …$\mspace{20mu} {{{\angle \; {{GM}_{i}\lbrack\; \ldots \;\rbrack}} = {{acos}\;\lbrack \frac{{MF}_{i} \cdot {GF}_{i}}{{{MF}_{i}}*{{GF}_{i}}} \rbrack}},{i = 1},2,\ldots}$  end$\mspace{20mu} {{\min \; \angle \; {TV\_ GM}_{i}} = {\overset{\_}{\angle \; {GM}_{\iota}} - {{tol}\; 1}}}$$\mspace{20mu} {{\max \; \angle \; {TV\_ GM}_{i}} = {\overset{\_}{\angle \; {GM}_{\iota}} + {{tol}\; 2}}}$  Where:  tol 1  and  tol 2  are  tolerances  for  the  detection  threshold.

A magnetometer radii check calculation is performed on the calibratedmagnetometer data before alignment. The radii of a normal point setshould be equal to one. Algorithm (7) may be used to determine the radiithreshold. As part of using Algorithm (7), a set of artificial points ona unit sphere are created or generated, where the set of points shouldbe close to evenly distributed and have a minimum of 18 points. Using anellipsoid fitting method on the artificial points, one can utilize theprocessed magnetometer data to determine the radii as set forth below

Radii=ellipsoidFit([artificial points,procMAG _(i) ], i=1,2 . . .

minRadii=Radii−tol3

maxRadii=Radii+tol4  Algorithm (7)

Where: tol3 and tol4 are tolerances for the detection threshold; and,

minRadii and maxRadii should be within 0.93 and 1.07 for applicationsthat require high reliability.

Quaternion output discrepancies is a detection method that is based onthe orientation output calculated from using different magnetometers.The following algorithm, Algorithm (8), may be used to determine thequaternion output discrepancies threshold, using the quaternions outputfrom Algorithm (6), and calculate the angle between the quaternions.

qDiff_(i,j)[ . . . ]=quatAngleDiff(q _(i) [ . . . ],q _(j)[ . . . ]),i≠j, i=1,2 . . . , j=1,2 . . .

TV_qDiff_(i,j)= qDiff_(i,j)[ . . . ]+SSF7*std(qDiff_(i,j)[ . . . ]),i≠j, i=1,2 . . . , j=1,2  Algorithm (8)

After extracting the TVs for anomaly detection, the followingcalculation may be used to determine if the IMU is being affected by alocal anomaly. As a prefatory step, all incoming sensor data may beprocessed using Equations 1-10 discussed previously. In addition,Algorithms (1) and (4) may be utilized to determine relative deviationamong magnetometers within a preset kernel size. For the exemplarysensor fault detection system, the motion is detected based on thefollowing rules in Table (1).

TABLE (1) Motion Detection based on sensors Activity Condition ResultGyro has detected If std(procGyro_(i)) > 

 , i = inertial_motion_detection = 1 motion 1, 2, . . . OR |procGyro_(i)| > 

 , i = 1, 2, . . . Accelerometer has If std(procAcc_(i)) > TV_(ACC) _(i), i = 1, 2, . . . inertial_motion_detection = 1 detected motionMagnetometer has If std(procMAG_(i)) > TV_(MAG) _(i) , i = 1, 2, . . .Magnetic_motion_detection = 1 detected motion

Within the preset kernel size, the following rules in Table (2) are usedto detect a motion anomaly.

TABLE (2) Motion anomaly detection Motion Anomaly Flag Condition FalseIf Σinertial_motion_detection > 1 AND Σmagnetic_motion_detection > 1False If Σinertial_motion_detection < 1 AND Σmagnetic_motion_detection <1 True If Σinertial_motion_detection < 1 ANDΣmagnetic_motion_detection > 1 False If Σinertial_motion_detection > 1AND Σmagnetic_motion_detection < 1

The following rules in Table (3) are used to detect a magnetic anomaly.

TABLE (3) Magnetic Anomaly detection Magnetic Anomaly Flag ConditionTrue costFunction > TV_(costFunction) False costFunction ≦TV_(costFunction)

In accordance with the instant disclosure, one can calculate themagnetic strength using Algorithm (5). Thereafter, using the rules ofTable (4), one can discern whether a magnetic anomaly is present.

TABLE (4) Magnetic strength anomaly detection Magnetic Strength AnomalyFlag Condition True MagStr_(i) < min_MagStr_(i) OR MagStr_(i) >max_MagStr_(i) False MagStr_(i) ≧ min_MagStr_(i) OR MagStr_(i) ≦max_MagStr_(i)

In accordance with the instant disclosure, one can determine the anglebetween magnetometer and accelerometer vectors using Algorithm (6). Thisdetection is only activated when the IMU/sensors is/are stationary,which is detected from the motion detection algorithm. Using thefollowing table, Table (5), one can discern whether a magnetometer andaccelerometer vector anomaly is present.

TABLE (5) magnetometer and accelerometer vectors anomaly detection MAGand ACC vectors Anomaly Flag Condition True ∠GM_(i) < min∠TV_GM_(i) OR∠GM_(i) > max∠TV_GM_(i) False ∠GM_(i) ≧ min∠TV_GM_(i) OR ∠GM_(i) ≦max∠TV_GM_(i)

In accordance with the instant disclosure, one can determine themagnetic radii by processing the raw data from the magnetometers usingAlgorithm (7). Using the following table, Table (6), one can discernwhether a magnetic radii anomaly is present.

TABLE (6) Magnetic Radii anomaly detection Magnetic Radii Anomaly FlagCondition True radii < minRadii OR radii > maxRadii False radii ≧minRadii OR radii ≦ maxRadii

In accordance with the instant disclosure, one can calculate theorientation of the IMUs with Algorithm (6) and calculate the angulardifference with Algorithm (8). Using the following table, Table (7), onecan discern whether a quaternion output discrepancy is present.

TABLE (7) Quaternion discrepancies detection Quaternion DiscrepanciesFlag Condition True qDiff_(i,j) > TV_qDiff_(i,j) False qDiff_(i,j) ≦TV_qDiff_(i,j)

Anomaly detection can be weighted using the output from each of thesesub-detection systems, and determine the relevant level of disturbancedetected by the system. For systems requiring high levels ofreliability, it should be considered an anomaly if any of thesesubsystem flags is detected.

The following discussion describes a process for calibrating an IMU whenit is being used with or in proximity to a ferromagnetic object. Oneexample, in the medical device field of use, is during surgicalnavigation, where the IMU is needed to track the motion of aferromagnetic object (e.g., tools made from stainless steel, CoCr, etc).Since local environment calibration is performed presumably withoutthese ferromagnetic objects in proximity, there is no way to compensatefor their distortions during the local environment calibration stage. Asecondary calibration step should be undertaken to correct thedistortion when ferromagnetic objects will be used in proximity to theIMUs. The proposed method accounts for the fact that an IMU may berigidly fixed to the ferromagnetic object when in use.

While it is possible to manually maneuver the IMU attached to theferromagnetic object to collect a lump sum of data for re-calibration,this method is inefficient and requires a substantial amount of time forthe user to repeat the maneuver for each ferromagnetic object needingcalibration. Secondly, it is possible that the calibration of thecombined effect of the local environment and the ferromagnetic object isnot aligned with the calibration of the local environment alone. Thougheither of the calibrations may be functional for yaw tracking, thelatent transformation between calibrations on the reference field canproduce errors. By way of example, one may think of it as comparing thereading of two compasses, one being horizontal to the ground, the otherbeing tilted 5 degrees. The horizontal compass is calibrated to thelocal environment, while the other compass is calibrated with bothferromagnetic objects and local environment. While both calibrations areboth functional, the reading from each compass will be different.Accordingly, an exemplary algorithm of the instant disclosure addressesboth of these issues.

The instant approach makes use of the distortion characteristic of lowfield strength ferromagnetic objects, and models the correction based ona limited amount of data input. This enables a much simpler calibrationmaneuver, while still being able to achieve reasonable calibration.

Referring to FIGS. 1-3, the object 100 attached to the IMU 102 (or theIMU 102 by itself) is moved through a series of orientations—which canbe confirmed from the IMU data. From the collected data, an initialestimation of the distortion field can be calculated through standardellipsoid processing and analysis. A reference dataset (currently using25000 vectors) is drawn from a von-Mises Fisher density, therebycreating a set of uniformly distributed vectors around the IMU sphere.The magnetic distortion calculated during the local environmentcalibration, and the initial estimation of magnetic distortion with theferromagnetic object are applied in reverse to the reference dataset.This creates two datasets: (1) a reference dataset being distorted bythe local environment; and, (2) a reference dataset being distorted byboth the local environment and the attached ferromagnetic object 100 (orferromagnetic object in proximity to the IMU 102). Since all the pointsbetween these two datasets originate from the same reference dataset, apoint-to-point correspondence between local environment distortion, andlocal environment distortion with the ferromagnetic object isestablished. The correction for the ferromagnetic object is thenextracted from the latter using a recursive optimization algorithm withthe prior dataset.

This exemplary calibration process can be used for any number ofobjects, so that multiple objects which have latent magnetic fields canbe used with the magnetometer without compromising sensor accuracy, aswell as maintain a consistent reference field.

The following section outlines the algorithm for calibrating aferromagnetic object with a pre-calibrated IMU. This part assumes thatthe IMU has been calibrated according to the foregoing discussion. TheIMU, as part of this calibration sequence, should be rigidly attached tothe tool or otherwise have the tool in a constant position with respectto the IMU. The operator rotates the object for at least one revolutionalong each axis in a standard Cartesian system as depicted in FIGS. 1-3.Once the tool calibration dataset is collected, the following algorithmcan calculate and compensate the distortion.

As an initial matter, an initial guess is made as to the center offsetof the data set using Equation 13 below.

$\begin{matrix}{{SFe}_{{inst}_{i}},{{HFe}_{{inst}_{i}} = {EF}_{{fun}{({DATA}_{{mag}_{i}})}}},{i = 1},2,\ldots} & {{Equation}\mspace{14mu} 13}\end{matrix}$

Thereafter, a large number (e.g., 25,000) of random 3D unit vectors aregenerated and application of the inverse of the original and instrumentSFe on to them using Algorithm (9) below.

Iteration: for 1 . . . 25000

samples=rand3Dvectors

samples_(distorted) _(i) =samples*inv(SFe _(inst) _(i) ), i=1,2, . . .

samples_(normal) _(i) =samples*inv(SFe _(orig) _(i) ), i=1,2, . . .

End  Algorithm (9)

Using Equation 14, shown below, calculate the point correspondencetransformation.

T _(i) _(→N) =samples_(normal) _(i) \samples_(distorted) _(i) ,i=1,2,  Equation 14

Using Equation 15, shown below, calculate the point-to-point distortiondistance.

Diff_(ND) _(i) =samples_(normal) _(i) −samples_(distorted) _(i) ,i=1,2,  Equation 15

Using Equation 16, shown below, for each axis, locate the maximumnegative of Diff_(ND).

loc_(i)=minmax(Diff_(ND) _(i) ), i=1,2,  Equation 16

Using Equation 17, shown below, calculate the index location of loc.

The index location of loc is

ind_(i)=index(minmax(Diff_(ND) _(i) )), i=1,2,  Equation 17

Using Equation 18, shown below, calculate the center along the maximumdistortion distance axes.

c _(i)=min(loc_(i))+max(loc_(i)), i=1,2,  Equation 18

Using Equation 19, shown below, calculate the maximum correctionboundaries.

D _(bound) _(i) =samples_(distorted) _(i) ,(ind_(i)), i=1,2,  Equation19

Using Equation 20, shown below, after the distortion compensationparameters are calculated, center the instrument calibration dataset.

centerMAG _(i)=DATA_(mag) _(i) −HFe _(inst) _(i) , i=1,2,  Equation 20

The foregoing calibration correction can be applied in at least twoways. In a first fashion, the calibration correction is applied using apoint correspondence transformation, such as Equation 21 shown below.

compensatedMAG _(i)=centerMAG _(i)*inv(T _(i) _(D→N) ) i=1,2,  Equation21

Alternatively, in a second fashion, the calibration correction isapplied using geometric scaling and compensation, such as by usingEquations 22 and 23.

$\begin{matrix}{{{S\; {F_{i}( {X,Y,Z} )}} = \frac{\begin{matrix}{{{centerMAG}_{i}( {X,Y,Z} )} - {{c( {X,Y,Z} )} \cdot}} \\{{{Dbound}_{i}( {X,Y,Z} )} - {c( {X,Y,Z} )}}\end{matrix}}{\begin{matrix}{{{norm}( {{{centerMAG}_{i}( {X,Y,Z} )} - {c( {X,Y,Z} )}} )} \cdot} \\{{norm}( {{{Dbound}_{i}( {X,Y,Z} )} - {c( {X,Y,Z} )}} )}\end{matrix}}},\mspace{20mu} {i = 1},2,\ldots} & {{Equation}\mspace{14mu} 22} \\{{{ccompensatedMAG}_{i} = {{centerMAG}_{i} + {{loc}_{i}*{SF}_{i}}}},{i - 1},2,\ldots} & {{Equation}\mspace{14mu} 23}\end{matrix}$

Using equation 24, the adjusted calibration parameters may becalculated.

SFe _(adjInst) _(i) ,HFe _(adjInst) _(i) =magCal(compensatedMAG _(i)),i=1,2,  Equation 24

Post conclusion of the foregoing calibration computations, thecalibration computations may be applied to the magnetometers of each IMUthat is attached to the ferromagnetic object or will be in a fixedpositional relationship thereto. In order to apply the calibrationcomputations to the magnetometers of each IMU, the magnetometer data iscentered using Equation 20. By way of example, on may apply thecorrection to account for the distortion cause by the tool via a pointcorrespondence transformation using Equation 21, or apply the geometricscaling and compensation using Equations 22 and 23. Thereafter, themagnetometers of the IMU can be processed by applying the adjusted softand hard iron compensation parameters via Equation 25, shown below.

$\begin{matrix}{{{procMAG}_{i} = {{processMag}( {{compensatedMAG}_{i},{SFe}_{{adjInst}_{i}},{HFe}_{{adjIinst}_{i}}} )}},\mspace{20mu} {i = 1},2,\ldots} & {{Equation}\mspace{14mu} 25}\end{matrix}$

Once the magnetometers have been processed, the magnetometers arealigned using Equations 9 and 10

Following from the above description, it should be apparent to those ofordinary skill in the art that, while the methods and apparatuses hereindescribed constitute exemplary embodiments of the present invention, theinvention described herein is not limited to any precise embodiment andthat changes may be made to such embodiments without departing from thescope of the invention as defined by the claims. Additionally, it is tobe understood that the invention is defined by the claims and it is notintended that any limitations or elements describing the exemplaryembodiments set forth herein are to be incorporated into theinterpretation of any claim element unless such limitation or element isexplicitly stated. Likewise, it is to be understood that it is notnecessary to meet any or all of the identified advantages or objects ofthe invention disclosed herein in order to fall within the scope of anyclaims, since the invention is defined by the claims and since inherentand/or unforeseen advantages of the present invention may exist eventhough they may not have been explicitly discussed herein.

What is claimed is:
 1. A method of calibrating an inertial measurementunit, the method comprising: collecting data from the inertialmeasurement unit while stationary as a first step; collecting data fromthe inertial measurement unit while repositioning the inertialmeasurement unit around three orthogonal axes of the inertialmeasurement unit as a second step; calibrating a plurality of gyroscopesusing the data collected during the first step and the second step;calibrating a plurality of magnetometers using the data collected duringthe first step and the second step; calibrating a plurality ofaccelerometers using the data collected during the first step and thesecond step; where calibrating the plurality of magnetometers includesextracting parameters for distortion detection and using the extractedparameters to determine if magnetic distortion is present within a localfield of the inertial measurement unit.
 2. The method of claim 1,wherein the act of extracting parameters for distortion detection adoptsat least one of the following assumptions: (1) motion of the inertialmeasurement unit is detectable by all of the plurality ofaccelerometers, the plurality of gyroscopes, and the plurality ofmagnetometer comprising the inertial measurement unit; (2) an angle anda length of the plurality of magnetometers have approximately identicalvalues; (3) each of a plurality of magnetometer vectors is approximatelyunity in magnitude; (4) an angle between a calibrated magnetometervector and a calibrated accelerometer vector should not change; (5) aradius of the plurality of magnetometers post calibration isapproximately equal to one; and, (6) quaternions calculated usingdifferent ones of the plurality of magnetometers are approximatelyequal.
 3. The method of claim 2, wherein a deviation from the assumptionthat motion of the inertial measurement unit is detectable by all of theplurality of accelerometers, the plurality of gyroscopes, and theplurality of magnetometer comprising the inertial measurement unit,results in a determination of a motion anomaly.
 4. The method of claim3, wherein a deviation from any one of the following assumptions resultsin a determination of a magnetic anomaly: (2) an angle and a length ofthe plurality of magnetometers have approximately identical values; (3)each of a plurality of magnetometer vectors is approximately unity inmagnitude; (4) an angle between a calibrated magnetometer vector and acalibrated accelerometer vector should not change; (5) a radius of theplurality of magnetometers post calibration is approximately equal toone; and, (6) quaternions calculated using different ones of theplurality of magnetometers are approximately equal.
 5. The method ofclaim 1, wherein calibrating the plurality of gyroscopes includesdetermining a bias of each gyroscope by arithmetic mean.
 6. The methodof claim 1, wherein calibrating the plurality of magnetometers includesdetermining whether at least one of any hard iron distortion is presentand any hard iron distortion is present.
 7. The method of claim 6,wherein determining whether at least one of any hard iron distortion ispresent and any hard iron distortion is present includes utilizing amulti-dimensional ellipsoidal fitting function to the data collectedfrom the plurality of magnetometers.
 8. The method of claim 1, whereincalibrating the plurality of magnetometers includes accounting foralignment differences among the magnetometers from a referencealignment.
 9. The method of claim 1, wherein calibrating the pluralityof magnetometers includes calibrating the plurality of magnetometers toa new reference alignment.
 10. The method of claim 1, whereincalibrating the plurality of accelerometers includes applying amulti-dimensional ellipsoidal fitting function to the data collectedfrom the plurality of accelerometers.
 11. The method of claim 1, whereincalibrating the plurality of gyroscopes includes extracting parametersfor distortion detection and using the extracted parameters to determineif a motion anomaly is present.
 12. The method of claim 1, whereincalibrating the plurality of accelerometers includes verifying that thedata collected from the plurality of accelerometers in the second stephas a motion value greater than the data collected from the plurality ofaccelerometers in the first step.
 13. The method of claim 1, whereincalibrating the plurality of gyroscopes includes verifying that the datacollected from the plurality of gyroscopes in the second step has amotion value greater than the data collected from the plurality ofgyroscopes in the first step.
 14. The method of claim 1, whereincalibrating the plurality of magnetometers includes verifying that thedata collected from the plurality of magnetometers in the second stephas a motion value greater than the data collected from the plurality ofmagnetometers in the first step.
 15. The method of claim 1, whereincalibrating the plurality of magnetometers includes applying amulti-dimensional ellipsoidal fitting function to the data collectedfrom the plurality of magnetometers.
 16. The method of claim 1, whereincalibrating the plurality of magnetometers includes calculatingorientation output quaternions from the data collected from theplurality of magnetometers.
 17. The method of claim 16, whereincalibrating the plurality of magnetometers includes calculating an anglebetween the calculated orientation output quaternions.
 18. The method ofclaim 1, wherein calibrating the plurality of magnetometers includesapplying a magnetometer radii calculations to the data collected fromthe plurality of magnetometers to determine if a magnetic anomalyexists.
 19. The method of claim 18, wherein if the magnetic anomalyexists, calculating a strength of the magnetic anomaly.
 20. An inertialmeasurement unit calibrated according to claim
 1. 21. An inertialmeasurement unit incorporated as part of a surgical navigation systemcalibrated according to claim
 1. 22. (canceled)